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Abstract. Wave self- focusing in molecular systems subject to thermal effects, such as thin 
molecular films and long biomolecules, can be modeled by stochastic versions of the Discrete 
Self- Trapping equation of Eilbeck, Lomdahl and Scott, and this can be approximated by 
continuum limits in the form of stochastic nonlinear Schrodinger equations. 

Previous studies directed at the SNLS approximations have indicated that the self- 
focusing of wave energy to highly localized states can be inhibited by phase noise (modeling 
thermal effects) and can be restored by phase damping (modeling heat radiation). 

Here it is shown that the continuum limit is probably ill-posed in the presence of spatially 
uncorrelated noise, at least with little or no damping, so that discrete models need to be 
addressed directly. Also, as has been noted by other authors, omission of damping produces 
highly unphysical results. 

Numerical results are presented here for the first time for the discrete models including 
the highly nonlinear damping term, and new numerical methods are introduced for this 
purpose. 

Previous conjectures are in general confirmed, and the damping is shown to strongly 
stabilize the highly localized states of the discrete models. It appears that the previously 
noted inhibition of nonlinear wave phenomena by noise is an artifact of modeling that 
includes the effects of heat, but not of heat loss. 



Foreword. This is a revised version of the material presented in several talks, given at the 
FPU+50 conference in Rouen, June 2005, at the University of Arizona, in September and 
November 2005, and at the University of New Mexico in October 2005. 

1. Introduction 

Combinations of mildly nonlinear wave motion in molecular structures with localized excita- 
tion modes can lead to localization of wave energy. Perhaps the early mathematical example 
was Davydov's soliton theory modeling exciton waves in protein molecules interacting with 
localized phonons, vibrations at CO bonds, where a continuum limit gave the integrable one 
dimensional focusing cubic Nonlinear Schrodinger (NLS) equation [5, 7, 6]. This was fur- 
ther developed by various authors including Scott [18], and extended to vibrations in other 
molecular systems such as crystalline acetanilide and in smaller molecules such as benzene 
[10, 11]. Approximations that eliminate the fastest internal vibration modes again lead to 
systems that are discrete counterparts of the ID focusing cubic NLS, or coupled systems of 
such. 

Two dimensional molecular structures such as Scheibe aggregates [3] lead to similar math- 
ematical models related to the 2D focusing cubic NLS [14, 2] . 
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Stochastic terms are a natural refinement, modeling effects such as random spatial variations 
in the medium (fixed pattern noise: time independent) or thermal agitation (white noise: 
temporally uncorrelated) [1]. Nonlinear optics has also produced continuum or semi-discrete 
examples including intense CW lasers and multi-cored optical fibers with random imperfec- 
tions in the medium or in the strength of the coupling between signals in the different cores 
or different propagation modes. 

However, noise without balancing losses can lead to an unphysical degree of spatial disorder 
or thermal runaway, impeding wave propagation in contradiction to experimental observa- 
tions [4]. Even worse, the obvious continuum limits give Stochastic NLS (SNLS) equations 
which seem to be well-posed only when the noise has adequate spatial correlation [8] , which 
is not necessarily consistent with the length scales in the molecular systems. In other words, 
noise can destroy the spatial smoothness needed to justify a continuum limit. 

More realistic modeling thus requires damping losses, to make possible attainment of ther- 
mal equilibrium, solutions with enough spatial smoothness to sustain traveling waves, and 
perhaps to justify a continuum limit model. This leads [4] to the Damped Stochastic Discrete 
Nonlinear Schrodinger equation (DSDNLS) 

(1) i^f + E J ^m + |*n| 2CT *n + 7*n^p - A*„^(|*„(i) \**) = 0, 

m 

with a = 1 giving the cubic case typical in physical applications. This is the main object of 
this study, but possible continuum limits will also be discussed briefly. The most obvious is 
the Damped Stochastic Nonlinear Schrodinger equation (DSNLS) 

(2) if + |*r*H*§! -A^«=0, 

derived and studied in [4]. Here W n (t;u) and w(x, t;u) are noise processes with u labeling 
realizations, and the physical meanings of the variables and parameters will be explained in 
the next section. 

The next section surveys the background for this study: physical origins, mathematical mod- 
eling, and previous results, some theoretical, but mostly numerical simulations. Section 3 
decribes the numerical methods used, and Section 4 presents numerical results, primarily 
for the surrogate case of a one dimensional structure with quintic nonlinearity a = 2. Then 
continuum limits are discussed in Section 5, with some suggestions as to how to overcome 
problems with previous approaches, followed by conclusions and a discussion of directions 
for further work. 



2. Background 



We give here some details of a two dimensional example from chemistry and biochem- 
istry: modeling Scheibe aggregates, a class of highly ordered thin films of molecules coupled 
by dipole interactions predominantly within a plane: an essentially two dimensional wave 
medium. However much of the modeling is also relevant to a variety of other molecular 
systems such as the essentially one dimensional protein models discussed above; see [6, 18] 
and references therein. 



2.1. Scheibe aggregates: highly ordered thin molecular films. Scheibe aggregates 
are highly regular arrangements of molecules in thin films, sometimes a single molecule thick, 
or with only weak interaction between neighboring layers of molecules. These structures have 
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important biological roles such as in photo-chemical reactions, and one laboratory example 
is the cyanine dye Scheibe aggregates first studied by Biicher, Kuhn, Mobius et al [3, 15]. 
They establish an arrangement of the molecules in a single layer "brick-wall" lattice, with the 
dominant dipole interactions being those with six nearest neighboring molecules, arranged 
in a hexagonal array of approximately dihedral D2 symmetry (half turns and reflection in 
two perpendicular axes). 

The molecules also have internal excitation states, which are coupled to the excitons and 
are also subject to thermal effects: random external forces due to collisions with molecules 
from outside the thin film. 



2.2. An exciton-phonon system with phase noise and damping. Such thin films can 
be modeled as an exciton-phonon system with noise and damping acting on the internal 
modes as described by Bang, Christiansen, If, Rasmussen and Gaididei in [1], which adds 
damping to the purely quantum mechanical modeling of Bartnick and Tuszyhski [2]: 

(3) ^~lr + ^ Jnm ^ m + X u n^n = 

(4) M^ + MX^ + MtiW-^A^-,^ 
where 

ty n (t) is the exciton wave at location n, 

u n (t) is the elastic degree of freedom of the molecule at location n, 
J nm is the dipole-dipole interaction energy, 
X is the exciton-phonon coupling constant, 

dB n /dt is random external forces, the formal time derivative of an independent Wiener 

process (Brownian motion) at each node, 

7 is the strength of the random external forces, 

A is the damping coefficient, 

M is the molecular mass, and 

fio is the Einstein frequency of each oscillator. 

These equations conserves the energy J\f = Y^ n \^n\ 2 under the Stratonovic interpretation 
of the stochastic term as discussed below, at least on a fully infinite lattice or with suitable 
boundary conditions such as periodic. 

The form of Eq's (3), (4) also covers a wide range of other applications such as the one 
dimensional protein molecule models mentioned above: the indices will often be taken below 
to simply enumerate a collection of nodes, with details such as spatial relationships between 
locations encoded in the coupling terms J nm . For further mathematical flexibility, the 
coupling term will henceforth be written as a general power law x| x I , n | 2f7 , though all physical 
models we know of have the cubic nonlinearity a = 1. The underlying spatial dimension 
will be denoted by D, so D = 2 in the model above. 

Without noise and damping, the system in (3), (4) is Hamiltonian, giving a second conserved 
quantity 

H = - Jnm*n*rn- rj^^U n \* n \ 2a + ^J2(< 2 + n U n)- 

n,m^n n n 
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2.3. Eliminating the phonon terms u n . The phonon terms u n can be eliminated using 
the variation of parameters formula which gives an integral expression for u n in terms of *f> n 
[4]. To eliminate the resulting time integral and initial data transients, one must restrict 
to times Xt >> 1 and make the slowly varying envelope approximation: that the exciton 
intensity |\I/ n (t)| 2 is slowly varying relative to the phonon frequency Slo- Thus, the presence 
of damping (A > 0) is essential. 

The reduced system is 

(5) = ih^+ Jnm^m + V\^n\ 2a ^n 

7Y T dW n XV d rl T , ,,, lT 

*n-^ ~ 7W^n(t)\ 2 }V n , where 



Mm dt Ql dt 
V = x 2 /Mn%, Q 2 = Ql - (A/2) 2 , and 

(6) -p- = X / e- Xs/2 sm(n S )B n (t-s)ds. 

dt Jo 

Note that the new noise processes W n are temporally correlated, except in certain limits 
such as strong damping. 

2.4. Rephrasing as a Damped Stochastic Discrete Nonlinear Schrodinger Equa- 
tion. If the coupling J nm is homogeneous, in that the quantities 

Jnn • — ^ ^ Jnm 

have a common value Jo, then the ^ n have a common average phase evolution e tJ °/ h . This 
can be removed by adding J nn to each coupling sum, and with some rescalings includ- 
ing 7x/(M/i 2 f2) —> 7, A/(M7g) — > A, one gets the Damped Stochastic Discrete Nonlinear 
Schrodinger equation Eq. (1). 

Even if the J nn are not all equal, one can use their average value J as the phase shift, and 
absorb the differences J nn — J into the noise coefficients as fixed pattern noise. 

The energy J\f = \^ n \ 2 is still conserved, and without noise or damping the system still 
has a conserved Hamiltonian, 



7Y--V/ * * Lrifl, 12(1+0-) 



1 + CT 

Simple examples are uniform nearest neighbor interaction on a line [ID] or square lattice 
[2D] of spacing I with all the non-zero J nm having the same value, J /I 2 . The coupling 
term is then a J times the standard three point second derivative [ID] or five point discrete 
laplacian [2D], and the equation is a discretization of the Damped Stochastic Nonlinear 
Schrodinger equation (2). Also, the Hamiltonian for the case of no noise or damping is the 
natural discretization using simple forward difference quotients for the gradient terms of the 
Hamiltonian for NLS, 



It can be useful in places to think of the ODE systems in relation to the familiar NLS 
equation, but a more careful consideration of continuum limit approximations is needed, as 
discussed in section 5. 
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2.5. Removing the extra time derivative term, and Stratonovic differential form. 

For some purposes, the time derivative should be eliminated from the damping term. Also, 
the rigorous mathematical formulation must be in terms of stochastic integrals and differ- 
entials, and in order to conserve the energy J\f, products involving stochastic terms must be 
interpreted in the Stratonovic sense. (Loosely, the Stratonovic integral is defined as the limit 
of midpoint rule (or trapezoid rule) approximations, whereas the ltd integral is the limit of 
left-hand end point Riemann sums: see [16] for details.) Solving for d*$> n /dt, substituting 
into the other time derivative leads to the stochastic differential form 



m \ m / 



(7) id^! n + 
with od denoting the Stratonovic differential. 



dt + -f^ n odW n = 



Why not convert to ltd integral form? Any system of stochastic ODE's in terms of the 
Stratonovic integral can be replaced by an equivalent system which gives the same solution 
under the ltd interpretation, by replacing the above Stratonovic differential by the corre- 
sponding Ito differential plus a correction term [16], and the Ito form is far more amenable 
to analysis such as existence and uniqueness proofs In the current case, this gives 

A 2 

(8) d^ n = [as before] + ijV n dW n - —^ n 

However, this form is undesirable for current purposes, particularly numerical simulations. 
The new term adds rapid exponential decay, destroying the manifestly conservative form. 
This reflects the fact that the Ito differential term itself generates rapid exponential growth 
of individual realizations, related to the fact that the ensemble average of Ito solutions 
satisfy the underlying noise-free equation, and so conserves the energy M. 

This prevents the use of time discretizations which inherently conserve energy M; such 
conservative discretizations are used here for the Stratonovic form. 

Also, the ltd form has no continuum limit with spatially uncorrelated noise. This might 
reflect the conjectured lack of existence of solutions to such continuum limits, even when 
the limit formally exists for the Stratonovic form. 

2.6. Wave self-focusing and energy localization in SNLS. In the continuum model 
of 2D NLS, self-focusing can lead to the formation of single point singularities, sometimes 
called wave collapse. The proof of this for the NLS is based on a variance argument, which 
has been extended to various cases of Stochastic NLS by Gaididei and Christiansen [13], 
Debussche and Di Menza [9], and Fannjiang [12], though all require noise that is sufficiently 
correlated in space and uncorrelated in time. The last mentioned author's results are as 
follows. 

For data of sufficiently rapid decay at infinity, the pulse width can be measured by its spatial 
variance 



V(t) = J |Vf ||x|| 2 dx. 



In the critical case of 2D cubic (and more generally, crD = 2), the ensemble average (V) is 
related to the ensemble average (Tt) of the Hamiltonian by 
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and in turn noise modifies conservation of the Hamiltonian to 

(10) ^ = /2:=i|*(p)||p|| a dp, 
where <$(p) is the power spectral density of the noise distribution. 

For spatially uncorrelated noise, $(p) is a positive constant, so R = oo: the formulas break 
down, but strongly suggest that SNLS has no solutions, at least in Sobolev space H . For 
comparison to the discrete equations, for which no analogues of these formulas are known, 
note that a discretization of SNLS with grid spacing I effectively has noise of correlation 
length scale Z, and such correlation in SNLS gives 

(11) R = o{J w )- 

Without noise, R = 0, (V) = V evolves quadratically, and TL < is sufficient condition for 
finite time singularity formation, as otherwise, V would become negative. Noise changes the 
evolution to the cubic 

(12) (V) = (F)(0) + bt + A(H)(0)t 2 + 2R/3t 3 . 

Clearly, with weak enough noise, Ti(0) < still leads to a prediction of negative (V) by some 
positive time to, so with positive probability, solutions must cease to exist before that time. 
However, sufficiently large noise eliminates this necessity, and hints at global existence, and 
at dispersion with a cubic rate of variance growth. 

Such formulas have not yet been extended to account for damping, but as seen below, there 
are hints of an additional negative term in dTt/dt. 

2.7. Wave self-focusing and energy localization in the discrete systems. With no 
noise or damping, numerical simulations of discrete counterparts of the NLS equation show 
phenomena analogous to wave collapse, even in the subcritical case of ID cubic, where the 
PDE has some degree of self-focusing, but cannot develop singularities. That is, solutions 
can have the energy concentrate until it is mostly at a single node (molecule), and then 
stay localized in a solution that seems to oscillate around a stable steady state. This was 
first described by Davydov [6] in models of protein molecules, and further analyzed and 
simulated by Eilbeck, Lomdahl and Scott [18, 10, 11], who describe the phenomenon as 
self -trapping. 

For the discrete systems with noise, no formulas are known for the evolution of ensem- 
ble averages considered above, but from Eq. (1) and alternative form (7) the Hamiltonian 
evolution for individual realizations satisfies 




These no longer make it clear that noise causes the previously noted linear increase in 
7i, or corresponding growth in beam spatial variance or inhibition of wave collapse, but all 
these are still seen in numerical studies of discrete systems with spatially uncorrelated noise, 
including those below. This is to be expected, since those discrete systems are effectively 
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discretizations of SNLS with noise of spatial correlation on length scale comparable to the 
mesh spacing of the discretization. 

The full 2D Stochastic NLS equation with spatially uncorrelated noise has been simulated 
by Bang, Christiansen, If, Rasmussen and Gaididei in[l], and the ID quintic case of the 
Stochastic NLS equation by DeBussche and Di Menza in [9]. In each case it is observed that 
spatially uncorrelated noise above a certain threshold level prevents wave collapse. 

The former paper indicates that this noise effect is too strong to match physical experiments: 
noise levels so low as to correspond to temperatures of a few Kelvin are needed to reproduce 
behavior seen in experiments at far higher temperatures. The likely cause is "thermal 
runaway" due to the absence of a mechanism for "heat loss" , such as the damping term. 

The latter authors also interpret simulations as showing that with spatially correlated noise 
(as is effectively imposed by a fixed spatial discretization), collapse can only be delayed, but 
will always occur. They also offer a non-rigorous argument for this second conclusion. 

However we come to a different conclusion below. In [9], the authors judge that collapse 
has occurred whenever the maximum amplitude has grown by a factor of three at any one 
time, but this seems unreliable when noise is present, since then amplitude spikes this high 
at a single node can occur transiently, rather than as part of ongoing focusing events. A 
more reliable criterion for numerical detection of wave collapse and energy localization is 
that the energy at a single node exceeds some substantial fraction of all energy, persistently 
for a significant interval of time. In the simulations below, localization is manifested with 
a majority of all energy staying at one node until the end of the computer time, and thus 
persisting for at least as long as the rise time of the localization. 

As to damping effects, Eq. (13) indicates that damping has the opposite effect of noise, 
causing reduction of 7i. Combined with the expectation that (9) still holds approximately 
for discretizations of SNLS, this suggests the possible return of wave collapse, in the discrete 
form of concentration of energy near a single node. 

Christiansen et al [4] have done the only simulations known to the current authors of the 
model with damping of Eq. (2). They do this with further approximate by a small system of 
ODE's, first imposing radially symmetry and then using the method of collective coordinates. 
Such modeling has lead to some analytical results on self-focusing in the NLS equation, but 
with the stochastic terms, it still requires study primarily by simulation. 

They start with simulations without damping, corroborating the above described obser- 
vations about inhibition of wave collapse by sufficiently strong noise, leading instead to 
dispersion of initially concentrated wave energy. With damping added, they observe that 
the effect of noise can be reversed, leading to collapse where with noise alone it would not 
occur, as suggested by Eq. (13). Again there appears to be a threshold damping level for 
this to occur with given initial data and noise level. 

With this survey done, we are ready to consider new numerical methods and results of 
simulations based on Eq. (1). 



3. Numerical methods: a variant of Chang and Xu's iterative trapezoid 

With noise but no damping, and with homogeneous coupling and periodic boundary condi- 
tions, a Fourier split-step method could be used, as was done by Bang et al [1] for SNLS. 
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Instead, an implicit time discretization based on fixed point iterative solution of the trapezoid 
rule is used, similar to one described and analyzed by Chang and Xu [17]. It has several 
virtues: 

• it satisfies the needed Stratonovic interpretation (no need for the Ito correction 
term), 

• it conserves the exciton energy J\f, and 

• with no noise or damping, it conserves the hamiltonian 7i. 

The main disadvantage is the need for iterative solution (so that conservation is no longer 
exact, but still highly accurate), and the fact that it is difficult to go beyond simple fixed 
point iteration due to the existence of coupled nonlinear terms, leading to the time step size 
restrictions typical of an explicit method rather than the underlying implicit method and a 
large number of iterations needed when noise effects are strong. Since the coupling is strong 
relative to the other terms (as in the discretized NLS equations that these essentially are) 
the system is rather stiff, leading to the need for rather small time steps. 

However, arguably the short spatial scales of the noise mean that the time step size limi- 
tations of simple iterative solution method are also the natural time scales of the smallest 
relevant features, so there might be little room for time step increases without inaccurate 
handling of noise effects. 

That is, the stiffness of the ODE systems is probably an essential time scale that must be 
preserved in the modeling, including the fully discrete model used for numerical solution. 



3.1. Trapezoid method: conservative time discretization. Writing VE^ for the ap- 
proximations of ^ n (P), 5t = P +1 — P , and 5Wi/5t for the approximation of cr n (t), constant 
on P < t < P + , the scheme used for the cubic case is 



(14) 







+i 



5t 



+ 



6w> . l*£| 2 + |tt£ +1 | 2 



st 



+ 



+AIm tf£ J mn^m. + *n +1 £ J, 



mn A rn 



3+1 



For the case studied so far of totally uncorrelated noise, the noise components W& are 
independent with normal distribution 5Wl ~ N (0, 5t/l) , so that 5W£/5t ~ N (0, l/(l5t)). 
The scaling with I and nearest coupling strength J nm = I /I 2 are used for consistency with 
interpretation as the discretization of the DSNLS, Eq. (2). 

More generally, exact conservation of the Hamiltonian is achieved by discretizing the non- 



linear term | tyj \ p 1 \Pj as 



p+l 



3+1 12 



l^nl P+1 * 
X — 



3+1 



3 12 



2 



Thus for quintic nonlinearity, the form for this term is not the familiar trapezoid approxi- 
mation, but 



|^ n | 4 + l^lWl 2 + |^n +1 | 4 &n + 



,3+1 



WAVE ENERGY LOCALIZATION BY SELF-FOCUSING IN LARGE MOLECULAR STRUCTURES 9 



3.2. Trapezoid method: iterative scheme. The nonlinear implicit scheme above is 
solved using a simple fixed point iteration, eliminating implicit form even for the linear 
dipole coupling (discrete laplacian) terms, and thus avoiding simultaneous linear equation 
solving, contrary to the original Chang- Xu algorithm. The reason is the nonlinearity and 
stiffness of the damping term and the formally unbounded noise terms, which lead in prac- 
tice to similar "explicit scheme" time step size restriction St = 0(l 2 ) even with implicit 
handling of the linear coupling term. 

Writing for the k-th iterate, the initial approximation used for the new time step 
is '° = *nj; an d and subsequent iterates are given by solving the uncoupled linear 
equations 



(15) 



St 



+ 



i*ii 2 +i*i +1 ' fc|2 



+ 7" 



swi_ 

St 



+ Aim & n ]T J mn W m + ]T J mnH 



i+l,k 



*L + *i +1 ' fc+1 



3.3. Time step size. The fixed point scheme has a worst case convergence rate of K = nSt 
with k = max n ^ n | J nm |/2, so k = 2 /I 2 for the three point discrete second derivative 
used in the ID quintic case, giving a convergence condition St < 1/re, = l 2 /2. To minimize 
computational cost in the sense of minimizing expected iterations per unit time, the optimal 
choice of K balances 0{\/K) time steps per unit time and the 0(lne/lnK) iterations per 
time step needed to meet a given error tolerance e. 

Minimization of | \ne/(K lnK)\) gives K = 1/e, which is in fact observed to be optimal with 
spatially uncorrelated noise and no damping, a case where it will be seen below that "thermal 
runaway" puts significant amount of signal in the shortest length scales, realizing the worst 
case convergence rate. On the other hand, without noise or with damping, solutions are 
smoother and larger values K ~ 0.5 are most efficient, with K < 1 always necessary for 
convergence. 



4. Numerical results 

For computational efficiency, simulations have been done with a single computational space 
dimension, using two different approaches to this reduction of dimension. 

The main studies are done for the one dimensional quintic case D = 1, a = 2, with this 
somewhat unnatural nonlinearity power used for the sake of remaining in the critical case 
for collapse in NLS. Homogeneous nearest neighbor coupling is used, corresponding to dis- 
cretizing NLS with the the standard three point discretization of the second derivative. 
Nearest neighbor coupling makes sense in the physical models, since dipole interactions are 
very short range. It also makes little sense to use higher order spatial discretizations of 
Stochastic NLS equations, because the noise eliminates the higher order smoothness needed 
to make such discretizations more accurate. This case, without damping, was also studied 
by DeBussche and Di Menza [9], as discussed in section 2.7 above. 
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The second reduction used is imposing radial symmetry on the two dimension cubic NLS 
(D = 2, a = 1), and then again using standard three point discretization of spatial deriva- 
tives. This allows comparison to the results of Bang, Christiansen et al [1, 4] also discussed 
above. 

4.1. ID lattice with quintic nonlinearity and homogeneous nearest neighbor cou- 
pling. The initial data used in this section is always discretization of ^f n (0) = i^o(x n ) = 
1.1(1 + cosx n )/2 on a uniform periodic grid of n no d es equally spaced nodes in the period 
cell [— n, it]. This is chosen to be close to the "Townes soliton" R$ central to theory of NLS 
self-focusing, giving Hamiltonian TL just slightly negative so as to ensure self-focusing in 
NLS, while having energy N = Nipo just slightly above the minimum value J\f(Ro) needed 
for self-focusing to be possible (c.f. [19]). 

The default parameters are n no d es = 100, and noise strength 7 = 0.04, damping strength 
A = 0.002 when noise or damping are present at all, with the choice of the latter two values 
explained below. Also, results for a single "standard noise realization" are presented in 
many graphs, with corroboration by data from multiple realizations where appropriate. 

In the figures, a systematic curve color coding is used. Black is used only for initial data 
with functions of node index, and for averages over multiple realizations with functions of 
time. The color sequence blue, green, red, cyan, magenta, is used both for later times 

with functions of node index, and for successive realizations with functions of time. 

4.1.1. Self-focusing and energy localization without noise or damping. The time evolution 
of Eq. (1) in the ID quintic case without any noise or damping is illustrated in Fig's 1 and 
2, which show the distribution of energy |^ n (t)| 2 amongst nodes, for various times t. The 
first figure shows all 100 nodes, at three times before self-trapping occurs at t = t* ~ 2.6 
and two times afterwards. The second and all subsequent graphs of energy distribution are 
restricted to nodes near the self-trapping locus. 

Note that once the energy at any node passes about 10, and certainly once it passes 20, 
energy is largely concentrated on only a few nodes so that the solution is no longer accurate 
as a discrete approximation of NLS. Likewise, in all subsequent solutions with the current 
spacing of I = 50/tt, data past the time when the energy at any node first exceeds about 10 
should only be considered as accurate for the spatially discrete models. Indeed the maximum 
possible energy at one node is N '/I, which for the current initial data is approximately 45.3. 
Thus the values seen here of about 35 and above at a single node represent a clear majority 
of all energy concentrated at that node. 

The evolution of the degree of energy localization is shown in Fig. 3, as measured by the 
maximum energy at any one node as a function of time. Once energy localization has 
occurred, it is seen to persist, but with significant oscillations. The oscillations fit with the 
idea of the solution entering a neighborhood of an orbitally stable stationary state that is 
a center: the existence of such localized stable center stationary states has been proven in 
the minimal case of n no des = 2 by Eilbeck, Lomdahl and Scott [11]. 

Note that this oscillatory behavior is purely a property of the discrete system, as it only 
sets in after the maximum single node energy becomes too high for the numerical solutions 
to be relevant to the related PDE models. 

4.1.2. Inhibition of self-focusing by sufficient noise, without damping. The results here are 
much as seen in the simulations by various previous authors discussed in section 2.7. With 
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1 D quintic, 100 nodes 

40 i 1 1 1 1 1 1 



35 - 



30 - 



>, 25 - 
<D 




10 20 30 40 50 60 70 80 90 100 

node 

Figure 1. Energy distribution for ID quintic, no noise or damping. Times 
and curve colors t = 0,2, 2.5, 3, 3.5, 4 as noted in the text. 

1 D quintic, 100 nodes 

40 i 1 1 1 1 1 1 , 1 1 1 




40 42 44 46 48 50 52 54 56 58 60 



node 

Figure 2. As in Fig. 1 but restricted to nodes near the self-trapping locus, 
showing energy persistently concentrated almost entirely at a single node. 
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1 D quintic, 100 nodes 
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t 

Figure 3. Evolution of maximum single node energy for ID quintic, no 
noise or damping. 

low levels of noise, up to 7 ~ 0.03, focusing of energy to a single node followed by persistent 
localization still occurs. Higher noise levels of 7 > 0.04 inhibit self- focusing and energy 
localization, as shown in Fig. 4 for the single standard noise realization, and confirmed by 
multiple realizations. 

Another notable feature is the loss of the spatial smoothness that would be needed for 
continuum limit PDE modeling. It seems likely that this spatial disorder has the effect of 
inhibiting exciton wave propagation, and that this is the mechanism which prevents energy 
localization, by preventing the needed energy flux. 

The evolution of the degree of localization, or lack thereof, is shown in Fig. 5. Note that 
transient spikes to values more than four times the initial value occur, above all close to 
the focusing time t* ~ 2.6 noted above, but these are not indications of focusing towards 
localization. Continuing this solution for considerably more time never again reaches the 
maximum of about 4.5 seen at t ~ 3. 

4.1.3. Effects of adding both damping and noise. The effect of adding damping at various 
strengths A = 0.001,0.002,0.01,0.1,0.26,0.27 while maintaining the noise level of 7 = 0.04 
is summarized in Fig. 6 in terms of the maximum single node energy, and spatial structure 
is shown for the case A = 0.002 in Fig. 7. 

Small damping values (A < 0.001) cause a transient increase in collapse, but not enough to 
produce localization; instead, one eventually gets dispersion, as with damped noise. With 
damping above some threshold near A = 0.002, collapse proceeds to persistent localization 
of most energy at a single node. With the larger damping values A = 0.01, 0.1, the maximum 
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1D quintic, 100 nodes, noise=0.04 
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Figure 4. Energy distribution for ID quintic, noise 7 = 0.04, no damping. 
Color-time labeling as in Fig. 1. 



1D quintic, 100 nodes, noise=0.04 




01 23456789 10 

t 



Figure 5. Maximum single node energy for ID quintic, noise 7 = 0.04, no damping. 
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Figure 6. Maximum single node energy for ID quintic, noise 7 = 0.04, 
damping from A = 0.001 to 0.27. 

1D quintic, 100 nodes, noise=0.04, damping=0.002 
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Figure 7. Energy distribution for ID quintic, noise 7 = 0.04, damping A = 0.002. 
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single node energy quickly becomes almost constant at a value very close to 45.3, which as 
noted above means that almost all energy is at one node. 

Apparently solutions settle into a very close approximation of a steady state that has en- 
ergy almost completely localized. This is also true at the lower damping values for which 
localization is seen, but with far slower onset. For A = 0.002 maximum single node energy 
continues the rising trend seen, reaching 41 by t = 50 and 43 by t = 100. 

This strong spatial localization even for A = 0.002 is shown in Fig. 8, which gives data for 
five times after onset of localization. This and the previous graph show a form far closer 
to a steady state than for the undamped, noiseless case above. The strong oscillations seen 
previously are absent here, replaced only by far smaller fluctuations, as are inevitably caused 
by the noise. 

1D quintic, 100 nodes, noise=0.04, damping=0.002, after localization 

40 | , 1 1 1 1 1 , 1 , 1 
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FIGURE 8. As above but for times t = 3.4,3.8,4.2,4.6,5 after the onset of 
energy localization, showing the persistence of a near steady state. 

Fig. 9 shows the approach to a nearly stationary state over a longer time interval as indicated 
by maximum single node energy. It also gives the solution with the same damping but no 
noise: there is relatively little difference, indicating that damping is the dominant mechanism 
driving the solution towards a steady state, and that this mechanism is robust enough to 
be little perturbed by noise. 

A final observation is that as damping strength A is increased, the localization initially occurs 
earlier, then later, and finally, localization completely fails with a sharp transition between 
A = 0.26 and 0.27. Sufficiently strong damping apparently causes spatial smoothing which 
not only counteracts noise effects as before but also inhibits self-focusing. 
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1D quintic, 100 nodes, with damping=0.002, with and without noise=0.04 
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Figure 9. Maximum single node energy for ID quintic, damping A = 0.002, 
for both noise A = 0.04 and no noise. 

4.1.4. Connections to evolution of the Hamiltonian. There are indications in Eq's (9-13) 
that the evolution of the Hamiltonian TL could be related to the occurrence of self-focusing 
and localization of energy, as it is for the NLS, so this will be examined. 

First, it can be shown that the result of Eq's (10,11) apply at least qualitatively to the 
discrete system with noise but no damping, by considering the latter as a discretization of 
the Stochastic NLS with noise having correlation length scale proportional to the lattice 
spacing I. Since the precise constant of proportionality is not known, this will be done by 
checking first that (TL) grows roughly linearly in time, and then by observing that this linear 
growth rate is roughly proportional to re^ orfes , as suggested by Eq. (11). 

The evolution of (TL) is approximated by the black curve in Fig. 10. This is the average 
of results for four noise realizations shown by the colored curves, of which the blue curve 
corresponds to the single realization used in earlier graphs. It is seen that there is indeed 
roughly linearly growth in time, at an average rate of about 15. 

Increasing n no( i es with initial data and noise adjusted as for grid refinement in discretization 
of the SNLS, the growth rate of the Hamiltonian is seen in Fig. 11 to go from about 15 for 
nnodes = 100 to 100-120 for n nodes = 200 and 740-1000 for n nodes = 400. Though only a 
single realization is shown in each case, the growth rates are again approximately linear, 
and the changes in the rates are consistent with the predicted n^ odes scaling. As noted in 
section 2.6, this indicates the ill-posedness of SNLS with spatially uncorrelated noise. 

The addition of damping is predicted by Eq. (13) to at least partially offset the growth of 
TL caused by noise. It is natural to ask under what circumstances this effect is sufficient to 
bring TL back to negative values, and how the restoration of negative TL is related to the 
restoration of energy localization, and the answers appear to be positive. 
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1D quintic, 100 nodes, noise=0.04, four realizations 
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Figure 10. Hamiltonian for ID quintic, noise 7 = 0.04, no damping. Col- 
ored curves are for different noise realizations, the black curve is their mean. 

1 D quintic, 200 and 400 nodes, noise=0.04, no damping 
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Figure 11. Hamiltonian for ID quintic, noise 7 = 0.04, no damping. Blue 
is with 200 nodes, red is with 400, one realization for each. 
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In the less interesting case of damping insufficient to restore energy localization, TL initially 
grows roughly linearly as without damping, but then levels out to fluctuation around a 
significantly positive value. 

With damping sufficient to restore energy localization (A = 0.002, 7 = 0.04), Fig. 12 shows 
for each of four realizations there is initial roughly linear growth of TL at about the same 
rate 15 seen above, but this is followed by slowing of the growth, and then a sudden drop to 
negative values. These negative values quickly become very large and remain so, as shown 
for the standard noise realization in Fig. 13. For each realization, the time of this drop is 
fairly close to the time at which localization occurs, and indeed is driven by the sixth power 
nonlinearity term in TL taking on a large negative value when the energy at any one node is 
a substantial proportion of the total. 

1D quintic, 100 nodes, noise=0.04, damping=0.002, four realizations 
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Figure 12. Hamiltonian for ID quintic, noise 7 = 0.04, damping A = 0.002. 
Same noise realizations and curve color-scheme as in Fig. 10. 

The longer time behavior of TL shown in Fig. 13 is continued decrease, correlated to the 
move closer to a localized steady state that is also indicated by Fig. 9. The values of both 
maximum single node energy and TL seem to settle slowly near asymptotic values, but with 
persistent noise driven fluctuations. However, whereas maximum single node energy gets 
quite close to the extreme value M/l ~ 45.3, the value of TL does not get so close to the 
minimum possible value of TL m in = 2(N — iV 3 /6)// 2 « —511, though both extrema occur for 
the same case of having all energy at one node. The extra energy is presumably persistent 
spatial variations or thermal energy in the free energy part of TL. 

4.2. Discrete 2D radially symmetric DSNLS. A nearest neighbor discretization of the 
2D cubic DSNLS with radial symmetry has also been studied, in order to test the conjectures 
described in section 2.7, based on a collective coordinates reduction of this 2D cubic case. 
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1D quintic, 100 nodes, noise=0.04, damping=0.002 
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Figure 13. The first realization above to t = 100. 

Another goal is computationally efficient initial comparison to simulations of the full 2D 
cubic model with noise but no damping by Bang et al [1]. 

A weakness is the physically unlikely assumption of radially symmetric noise, but the qualita- 
tive features match those seen in the previous section, suggesting robustness under variations 
of lattice geometry, coupling, and nonlinearity power. 

Indeed, no fundamental differences are seen from the previous study of D = 1, a = 2 case, 
so we will simply summarize the points of agreement. 

All results are for the initial data from discretization of ^(0, r) = 1.65sech(r) on < r < 10 
with 256 equally spaced nodes. This function is close to the ground state "Townes soliton", 
giving initial value of the hamiltonean Ti ~ —0.006 and J\f ~ 1.8867, just above the collapse 
threshold value N c ~ 1.87. Thus singularity formation is ensured in the absence of noise 
and damping, but only just. The idea is to maximize the sensitivity of self-focusing effects 
to the perturbations applied. 

Noise without damping inhibits self-focusing localization of energy if its strength 7 exceeds 
a threshold, between 0.01 and 0.02. The hamiltonian for each realization Ti grown roughly 
linearly in time, and more so the approximation of ensemble average {Ti) given by averaging 
even a few realizations. If the spatial discretization is refined to more nodes, the rate 
constant increases rapidly, suggesting the failure of the continuum limit and ill-posedness of 
this form of SNLS. 

With above-threshold noise 7 = 0.02, damping restores energy localization once its strength 
A in turn exceeds a threshold, between 0.05 and 0.1. This is related to first slower growth 
of Ti, and then its decrease to substantially negative values at about the same time as the 
localization. 
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4.3. ID lattice, quintic nonlinearity, spectral coupling. Finally, a different spatial 
discretization has been used at the opposite extreme to nearest neighbor: discretization of 
the damped stochastic NLS equation for D = 1, a = 2 with periodic boundary conditions 
as above, but with spectral discretization of derivatives. 

The only changes seen are in quantitative details, not the qualitative features described for 
both previous cases, so details are omitted. 

This coupling is unnatural for the underlying molecular systems due to a mixture of signs 
corresponding to alternation between attractive and repulsive dipole interactions. However 
if continuum limits are valid for suitable combinations of noise and damping, so that there 
is smoothness on a scale larger than the molecular spacing, spectral discretization could 
potentially reduce computational costs. 

5. Refining the continuum limit approximations 

As discussed above in sections 2.6, 4.1.2 and 4.1.4, the continuum limit from the undamped 
Stochastic Discrete NLS (Eq. 1 with A = 0) with spatially uncorrelated noise is probably 
ill-posed; a so-called "ultra-violet catastrophe" . Some brief comments on possible solutions 
are offered to indicate possible future research directions. 

One possibility is that the addition of damping restores well-posedness, and numerical evi- 
dence above suggests that this might be so. However, it is not clear that the damping term 
allows one to overcome the technical obstacles to establishing existence and uniqueness; 
instead its fully nonlinear term might only increase the technical difficulties. 

A second approach is retaining higher order terms in approximating the discrete coupling 
terms by Taylor series expansions, which preserves explicit dependence on a length scale 
parameter I. 

Assuming the Dihedral D2 symmetry of the brick-wall molecular film structure, one gets 

(16) ^2 Jnm ^ m = hfitpxx + jo,2^yy + ^[jAfl^xxxx + h^xxyy + joA^yyyy)] + 0{l A ) 

in 

With the physically natural assumption of non-trivial and attractive coupling, meaning that 
all J nm are non-negative and some are positive, all the new j a b coefficients are non-negative 
and all j a o,job are positive. Linear rescaling of the x and y variables can then give the form 

(17) Yl J nn^m = ^xx + ^yy + 0(l 2 ), 

Til 

and discarding terms that vanish as I — > lead to the standard DSNLS approximation in 
Eq. (2). 

If instead we retain the next most important terms, involving fourth order spatial derivatives, 
and assuming D2 symmetry, one can rescale to get 

dtp d(\ifj\ 2 ) 

(18) i— + Alp+l 2 [j4fi(^xxxx + Ipyyyy) + 3 2,2^ xxyy] + + l<T1p ~ A ^ 1p = 0. 

With nearest neighbor coupling on a rectangular lattice the cross derivative term vanishes 
(J 2, 2 = 0), but that does not fit the brick- wall symmetry observed for cyanine dye Scheibe 
aggregates. 

A third approach is to use Pade fitting to the coupling term instead of Taylor polynomials, 
giving a pseudo-differential equation. With nearest neighbor coupling on a rectangular 
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lattice one can get 



One potentially advantage is that the pseudo-differential operator is bounded, which explicit 
removes the risk of an "ultra-violet catastrophe" and so might facilitate analysis. 



• Noise without damping can inhibit wave collapse, related to an increase in the hamil- 
tonian. 

• There is a minimum threshold noise level needed to do this. 

• The continuum limit as a Stochastic NLS equation with spatially uncorrelated noise 
appears to be ill-posed. 

• Damping can restore collapse, again requiring a threshold strength. 

• Damping might also restore well-posedness with spatially uncorrelated noise. 

• Even more damping inhibits self-focusing and collapse. 

In further work, the model should be refined in various ways. 

• Simulations and analysis of full 2D models. 

• Consideration of the ID cubic case and related models of nonlinear waves in long 
bio-molecules, where energy localization still occurs in discrete NLS models but the 
NLS continuum limit does not have self-focusing collapse. 

• Analysis of well-posedness of the various PDE and pseudo-differential equation mod- 
els, including the effects of damping. 

• Direct analysis and study of the ODE system (lattice models) with various interac- 
tion forms such as longer range. 

• Simulation with time correlation in the noise, as seen in Eq. (6). 
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